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ABSTRACT 

We use 3D-PDR, a three-dimensional astrochemistry code for modeling photodissociation regions 
(PDRs), to post-process hydrodynamic simulations of turbulent, star-forming clouds. We focus on the 
transition from atomic to molecular gas, with specific attention to the formation and distribution of H, 
C+, C, H2 and CO. First, we demonstrate that the details of the cloud chemistry and our conclusions 
are insensitive to the simulation spatial resolution, to the resolution at the cloud edge, and to the ray 
angular resolution. We then investigate the effect of geometry and simulation parameters on chemical 
abundances and find weak dependence on cloud morphology as dictated by gravity and turbulent Mach 
number. For a uniform external radiation field, we find similar distributions to those derived using a 
one-dimensional PDR code. However, we demonstrate that a three-dimensional treatment is necessary 
for a spatially varying external field, and we caution against using one-dimensional treatments for 
non-symmetric problems. We compare our results with the work of Glover et al. (2010), who self- 
consistently followed the time evolution of molecule formation in hydrodynamic simulations using a 
reduced chemical network. In general, we find good agreement with this in situ approach for C and CO 
abundances. However, the temperature and H2 abundances are discrepant in the boundary regions 
(Av < 5), which is due to the different number of rays used by the two approaches. 
Subject headings: astrochemistry, hydrodynamics, molecular processes, turbulence, stars: formation, 
ISMimolecules 



1. INTRODUCTION 

In the local universe, stars appear to form exclusively 
in cold, dense clouds of predominately molecular gas 
(McKee & Ostriker 2007). Understanding the evolution 
of these molecular clouds (MCs) and the formation of 
stars within them is a fundamental problem in astro- 
physics that is hampered by distance, projection effects, 
and the high optical depth in these regions. Probing the 
mass and velocity distributions of the gas is further com- 
plicated by the fact that the most abundant molecule, 
H2, lacks a dipole moment. The next most abundant 
molecule, CO, which is commonly used to probe the cold 
molecular gas distribution in lieu of H2 , has a typical av- 
erage abundance of about one per 10^ H2 molecules in 
the Milky Way. In addition, the relationship between CO 
abundance and total gas mass is a complicated one that 
depends upon metallicity, the three-dimensional radia- 
tion field, the abundances of other molecules, and dust 
chemistry (Beh et al. 2006; Glover & Mac Low 2011; 
Shetty et al. 2011). Accurately modeling the formation of 
H2 and the relative abundances of homologous molecules 
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such as CO requires following complex chemical reaction 
networks that encompass hundreds of species and thou- 
sands of reactions. 

Traditionally, the computational expense of evolving 
large chemical networks limited astrochemical investi- 
gations to simple one-dimensional hydrodynamic mod- 
els (e.g., Bergin et al. 2004) or to post-processing (e.g., 
Levrier et al. 2012). However, in recent years "re- 
duced" chemical networks have been adopted to investi- 
gate chemistry concurrently with three-dimensional hy- 
drodynamics (Nelson & Langer 1997, 1999; Pavlovski 
et al. 2002, 2006; Glover & Mac Low 2007a,b; Glover 
et al. 2010). Such methods have the advantage of being 
able to follow the temperature evolution of the gas due 
to UV heating and atomic and molecular cooling, which 
in principle infiuences the gas dynamics since shock jump 
conditions depend upon the local temperature. Nonethe- 
less, the expense of following the molecular evolution in 
situ necessitates various simplifications, including neglect 
of dust physics and coarse treatment of the radiation 
field. 

Thus far, turbulent cloud calculations including sim- 
plified chemistry have also focused on larger cloud com- 
plexes and generally neglected the self-gravity of the 



gas (see Glover & Clark 2012 as an exception includ- 
ing gravity). Neglecting gravity obviates the need for 
considerable additional resolution which would other- 
wise be required to resolve collapsing gas (Truelove et al. 
1997). In addition, without forming embedded sources 
to provide additional radiation (e.g., Offner et al. 2009; 
Krumholz et al. 2007), heating depends only on the ex- 
ternal cloud environment, leading to simpler radiative 
conditions. The gas temperature range induced by a 
standard external interstellar radiation field is generally 
limited (<100K) and deviates from 10 K mainly at low 
Av. 

Despite such simplifications, the astrochemistry under 
investigation is rich and not well understood. For exam- 
ple, cloud boundary regions are especially interesting be- 
cause this is where gas transitions from being ionized and 
atomic to predominantly molecular. These low-Av tran- 
sitions areas are by definition PDRs, where FUV photons 
dominate the energy balance and gas chemistry. PDRs 
are ubiquitous in the interstellar medium and are the 
source of most of the infrared radiation in galaxies. The 
recent development of 3D-PDR (Bisbas et al. 2012, here- 
after B12), which is the first dedicated FDR code able 
to treat arbitrary three-dimensional density distribution, 
now allows the accurate study of these regions in more 
complex structures. 

We dedicate this paper to three main goals. First, 
we compare 3D and ID treatments of a complex PDR 
region in order to evaluate the impact of dimensional- 
ity on chemical results. Thus, we extend the work of 
B12, who demonstrated the importance of higher dimen- 
sional treatment in accurately modeling simple 3D prob- 
lems, to consider complex, turbulent gas distributions. 
Second, we use self-gravitating, hydrodynamic simula- 
tions of molecular clouds with different Mach numbers 
to evaluate the importance of underlying physical pa- 
rameters on chemical abundances and distributions. Fi- 
nally, we explore the differences between two astrochem- 
istry approaches by considering results obtained via post- 
processing using 3D-PDR and results obtained from a 
chemical network calculation preformed "in situ" (e.g.. 
Glover et al. 2010). 

The paper is organized as follows. In section 2 we 
describe the 3D-PDR methodology and our hydrodynamic 
numerical simulations. In section 3 we validate our choice 
of spatial resolution by presenting convergence studies 
of grid-sampling in the cloud interior and at the cloud 
boundaries. We present our results in section 4, including 
a comparison to Glover et al. (2010) and discussions of 
chemical dependence on domain dimensionality, external 
radiation field, and cloud physical parameters. Section 5 
contains a discussion of future work and conclusions. 



2. METHODS 



2.1. 



Hydrodynamic Simulations 

In this paper, we analyze snapshots of four different 
hydrodynamic simulations of turbulent molecular clouds. 
The simulation parameters are summarized in Table 1. 
Three of the simulations (Rm4, Rm6 and Rm9) are per- 
formed with the ORION adaptive mesh refinement (AMR) 
code (Truelove et al. 1998; Klein 1999). Since these sim- 
ulations have not been previously published, we describe 
our method in detail below. 



ORION employs a conservative second order Godunov 
scheme to solve the equations of compressible gas dy- 
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where p, P, v are the gas density, pressure, and velocity, 
respectively. Here, e is the total energy e = ^pv^ + -^j, 
where 7 is the ratio of specific heats. ORION solves the 
Poisson equation for the gravitational potential, (j): 
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where rUn and Xn are the mass and position of the nth 
star, respectively. 

We close these equations with an isothermal equations 
of state: 

P = p , (5) 

where ks is the Boltzmann constant, /ip = 2.33 is the 
mean mass per particle, mu is the hydrogen mass, and 
T = 10 K is the isothermal gas temperature. Authors 
sometimes adopt a barotropic equation of state (e.g., 
Offner et al. 2008), which sets a characteristic density 
above which the gas becomes optically thick and ceases to 
be isothermal. However, the density at which this occurs, 
pc ^ 10"^"^ g cm~^, as calculated using full radiative 
transfer (Masunaga et al. 1998), exceeds the maximum 
density at our maximum AMR resolution (~ 5 x 10~^^ 
g cm~^). Consequently, the isothermal approximation is 
appropriate here. Alternatively, we might solve for the 
radiation field using a fiux-limited diffusion (FLD) ap- 
proach and thus take into account heating from forming 
stars (Offner et al. 2009). This would be more numeri- 
cally expensive but more physically accurate in the dense 
star- forming gas. However, without some prescription 
for protostellar outflows the stellar heating in the calcu- 
lation would be an over-estimate (Hansen et al. 2012), 
and moreover, an FLD approach would not supply more 
accurate information about the temperatures of the low- 
extinction gas as 3D-PDR does. 

We insert finer AMR grids when the local density vio- 
lates the Truelove criterion (Truelove et al. 1997): 



p<pj = r 
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where Axi is the cell size on level / and we adopt a Jeans 
number of J = 0.125. A sink particle is inserted when the 
gas exceeds the Jeans density for J = 0.25 on the maxi- 
mum level (Krumholz et al. 2004). In this paper, we do 
not analyze the sink particle distribution and properties; 
these are the subject of Kirk et al. (in preparation). 

We initialize the simulations with uniform density and 
then perturb the gas for three crossing times using a 
random velocity field (e.g., Mac Low 1999). This field 
has a fiat power spectrum for wavenumbers k = 1..2, 
which corresponds to physical scales of L..L/2. We 
re-normalize the perturbations to maintain a constant 
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Fig. 1. — Total gas column density for ORION snapshots Rm6_1.0 (top left), Rm9_1.0 (top right), Rm6_0.0 (bottom left), Rm4_1.0 (bottom 
right). The snapshot times are Itff , 1%, 0%, and 1%, respectively. 



cloud velocity dispersion. In the fiducial simulation, 
Rm6, the Mach number is chosen to satisfy the observed 
linewidth-size relation (McKee & Ostriker 2007). Fol- 
lowing the driving initialization, the simulations achieve 
a well-mixed turbulent state and we turn on gravity, al- 
lowing collapse to proceed for a global free-fall time. 

The ORION simulations all have a 256^ base grid and 
four levels of AMR refinement. As summarized in Table 
1, these three calculation have a total gas mass of 600 
M0, domain size of 2 pc {Ax^ = 100 AU), and turbulent 
3D Mach numbers of 4.2, 6.6 and 8.9. For comparison, we 
also analyze Rm6 without gravity, i.e., at t = Otff, and 
at half a free-fall time. Figure 1 shows the integrated 
column density at one free-fall time for these runs. 

We include the third simulation, n300, in order to di- 
rectly compare our PDR methodology to that of Glover 
et al. (2010), henceforth GIO. The n300 simulation was 
performed by S. Glover with a modified version of ZEUS- 
MP, which tracks the abundances of 32 chemical species. 
The n300 calculation uses a fixed 256^ grid. Turbulence 



is generated using random velocity perturbations in a 
manner similar to that used for the ORION simulations. 
It does not include self-gravity but does solve the equa- 
tions of ideal magneto-hydrodynamics and begins with 
an initially uniform magnetic field of 6 /iG. 

Figure 2 shows the mass-weighted and volume- 
weighted density distributions and corresponding chem- 
ical regimes for each of the ORION snapshots. The den- 
sity distribution functions exhibit a characteristic log- 
normal shape as expected for supersonic turbulent gas 
(e.g., Padoan et al. 1997; Kritsuk et al. 2007). As 
self- gravity becomes important, the density distribution 
grows a high-density tail (Mac Low & Klessen 2004). 
The cells at the peak of the density distribution fall into 
the PDR regime for the simulation parameters we adopt. 
The vertical lines in the histogram indicate the division 
between ionized, PDR and molecular gas. 

2.2. 3D-PDR 

3D-PDR (Bisbas et al. 2012) is a three-dimensional 
time-dependent astrochemistry code for treating pho- 
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Fig. 2. Density distributions for runs Rm9_1.0_12 (black, solid), Rm6_1.0_12 (red, dot), Rm4_1.0_12 (purple, dash), Rm6_0.5_12 (blue, 
dot-dash), and Rm6_0.0_12 (green, dot-dot-dash). The gas state is characterized as ionized (left), PDR (middle) or molecular (right), where 
vertical lines indicate the state boundaries. 

relative heating rates.) The thermal balance is solved 
self-consistently with the chemistry to determine the gas 
temperature. Unless otherwise specified, we adopt to- 
tal Carbon and Oxygen abundances of xq = 10~^ and 
xq =3.16 X 10~^. Further details can be found in B12. 

For the purposes of this paper we consider as PDR 
any H-nucleus density within the region 200 < nn < 
lO^cm"^. Below nn = 200 cm~^ we consider it ionized, 
whereas above nn = 10^ cm~^ we consider it fully molec- 
ular, with constant gas temperature and abundances that 
are independent of the external radiation field. The lower 
density limit is somewhat arbitrary since the H to H2 
transition can occur down to lower densities depending 
on the temperature. We impose this cutoff on the PDR 
calculations since we assume that gas at lower densities 
represents the HII component of the medium, which can 
only be reliably modeled using a photoionization code 
(e.g., MOCASSIN Ercolano et al. (2003, 2005, 2008). 

In this paper, once the gas is fully molecular we do 
not solve for its properties with 3D-PDR. Instead, we 
adopt the limiting values of the temperature and abun- 
dances for a uniform density of nn = lO^cm"^, which 
correspond to 10 K and nco/^H = 10~^, wherein no 
atomic Carbon remains. This is a reasonable approxi- 
mation for these densities since this gas, by definition, is 
well shielded from the external radiation and is almost 
entirely molecular. 

The cosmic ray ionization rate per H2 molecule is taken 
to be C = 5 X 10~^^ s~^ . The dust temperature is constant 
and set to Tdust = 20 K. We use A/^ = 12 rays of healfix 
refinement (level ^ = 0) and we use ^crit = 0.5 (c^ 7r/6) rad 
for the search angle criterion. We neglect the contribu- 
tion of the diffusive component of the FUV field by in- 
voking the on-the-spot approximation (Osterbrock 1974). 
We consider we have obtained thermal balance either 
when the heating and cooling rates differ by derr < 0.5%, 
or when the difference in temperature between two con- 
secutive iterations is Tdiff < 0.01 K. Finally, we typically 
evolve the 3D-PDR simulation to final times from 5.7—100 
Myr at which point the chemistry is in equilibrium (e.g., 
Bayet et al. 2009). Table 2 summarizes all the runs we 
perform with 3D-PDR. 



todissociation regions (PDRs) of arbitrary density dis- 
tribution. The code is able to solve self-consistently the 
chemistry and the thermal balance within any three- 
dimensional cloud. It uses an escape probability ap- 
proximation (or Large Velocity Gradient - Sobolev 1960; 
Castor 1970; de Jong et al. 1975) to compute the cooling 
functions. To do this, 3D-PDR uses a ray tracing scheme 
in which the directions of the rays are controlled by the 
HEALFIX algorithm (Gorski et al. 2005). This ray trac- 
ing scheme creates a discrete set of evaluation points by 
projecting the elements of the cloud along each ray. It 
can thus evaluate the column densities, the attenuation 
of the far ultraviolet radiation into the PDR, and the 
propagation of the FIR/submm line emission out of the 
PDR. 

As a further development of the fully bench-marked 
one-dimensional ucl_fdr code (Beh et al. 2006), 3d- 
FDR adopts the same chemical model features. For the 
simulations presented in this paper, we use a chemical 
network which is a subset of the UMIST data base of 
reaction rates (Woodall et al. 2007). This "reduced" net- 
work consists of 320 reactions and 33 species (including 
electrons). However, 3D-FDR also includes heating due to 
photoionization and photodissociation reactions in addi- 
tion to the standard gas-phase chemistry. Self-shielding 
of H2 and CO against photodissociation is accounted for. 
Comprehensive treatment of various gas heating mech- 
anisms (i.e., photoelectric heating from dust grains and 
PAHs, collisional de-excitation of vibrationally excited 
II2 following FUV pumping, photoionization of neutral 
carbon, cosmic ray heating) and emission from major 
cooling lines ([CII], [CI], [01], CO) are calculated at each 
element. 3d- FDR also includes turbulent heating, which 
is proportional to v^jj^^^/L^ where ^'turb is the turbu- 
lent velocity and L is the integral scale. Here, we adopt 
constant values of L = 5 pc and ^turb = lkms~^. In 
practice, L should be set to the simulation domain size 
and 'L'TURB to the ID turbulent Mach number times the 
mean sound speed, however we find that the turbulent 
heating is small compared to photoelectric, cosmic-ray 
and chemical heating, which are the other main sources 
of heating. (See the Appendix for a discussion of the 





TABLE 1 

Simulation Properties 




Snapshot ID^ 


L(pc) 


M(Mq) M A^min-A^max^ 


t/ts 


Rm6_0.0 
Rm6_0.5 
Rm6_1.0 
Rm9_1.0 
Rm4_1.0 
n300 


2 
2 
2 
2 
2 
20 


600 6.6 1..2 
600 6.6 1..2 
600 6.6 1..2 
600 8.9 1..2 
600 4.2 1..2 
82800 12.5 1..2 


0.0 
0.5 
1.0 
1.0 
1.0 
0.0 



^Simulation output ID, box length, total initial gas mass, Mach 
number, and fraction of a global free-fall time with gravity, re- 
spectively. 

The wavenumber range of the random velocity perturbations. 

Although Rm4, Rm6 and Rm9 each have 4 levels of 
grid refinement with a minimum cell size of 100 AU, 
we consider only the 256^ base-grid data when post- 
processing. The refined cloud regions, by construction, 
contain high-density gas that is > 10^ cm~^. At these 
densities, 3D-PDR considers the gas to be fully molecular 
and adopts a constant gas temperature and abundances. 

2.3. ''One-Way^' Hydrodynamic- Chemical Coupling 

Our method can be considered a "one-way" code cou- 
pling, because 3D-PDR uses the density output of the 
hydrodynamic calculations to compute the chemical dis- 
tribution. A benefit of this approach is that it is compu- 
tationally efficient, and large networks of reactions may 
be considered that would otherwise be too time consum- 
ing to compute in combination with the hydrodynamics. 
In addition, the affects of different radiative conditions 
and metallicity may be studied using the same hydrody- 
namical simulation. The deficit to this approach is that 
the corresponding temperatures computed by 3D-PDR do 
not affect the subsequent hydrodynamic evolution. In 
a one-way coupling, consistency between the hydrody- 
namic quantities and chemistry is only achieved if the a 
priori simulated values are chosen to reffect the antici- 
pated post-processed values. Because 3D-PDR computes 
a wide distribution of temperatures, it is not possible to 
achieve consistency by adopting a single, constant tem- 
perature. For example, for Rm6_1.0_12 3d-pdr deter- 
mines a mass-weighted temperature of ~22 K, which is a 
factor of two above the fiducial 10 K simulation tempera- 
ture. However, because we adopt 10 K for the simulation, 
by construction the densest regions, i.e., the star-forming 
gas (n > a few 10^), their dynamics will be in fairly good 
agreement with the computed 3d-pdr temperatures. It 
is also worth noting that for a simulation with a ID rms 
velocity of 0.7 km s~^, gas temperatures would need to 
reach ~ 140 K in order to obtain dynamic parity with the 
turbulent gas pressure (assuming a stellar external radia- 
tion field). Since the 3d-pdr computed gas temperatures 
are generally much less than 140 K, the hydrodynamics 
would remain governed by turbulence and so only a small 
difference would be expected if the 3d-pdr temperatures 
were fed back into the simulation. 

In the simulation we also adopt a fixed value for the 
mean mass per particle, /ip, which implicitly assumes 
that the gas is entirely molecular. We will show later that 
the hydrogen is almost all in molecular form through- 
out the domain with the exception of a few cells at the 
domain edge. Since molecular hydrogen dominates the 



mass budget of the gas by several orders of magnitude 
this particle mass approximation is a good one for the 
simulations used in this study. 

A second discrepancy between the dynamics and the 
chemistry occurs because 3D-PDR assumes that the ra- 
diation field impinges on the gas at the box boundaries, 
while the hydrodynamics assume periodic boundary con- 
ditions, i.e., there is no edge. This incongruity is also part 
of the GIO approach, which adopts periodicity for the 
gas but not the radiation field. For any boundary con- 
vention, high-density gas will have high-extinction nearly 
independently of location with respect to the boundary. 
Since turbulent clouds are naturally porous and the dense 
gas has a low- volume filling fraction, we can expect that 
radiation would penetrate many lower density regions for 
some sight-line to the "edge." Practically, the effect of 
the incident radiation field is to define a new effective 
boundary for the molecular gas, which reflects the flla- 
mentary and inhomogenous shape of the gas. Authors 
that seek to model an entire cloud rather than a periodic 
piece must instead wrestle with the arguably equally difli- 
cult problem of how the cloud connects to the larger-scale 
ISM, which is related to the issue of molecular cloud for- 
mation (e.g., Banerjee et al. 2009; Van Loo et al. 2013). 

3. METHOD VALIDATION 
3.1. Crid Sampling 

We flrst verify that our results are converged and in- 
dependent of the 3D-PDR grid resolution by compar- 
ing the calculated abundances for the same simulation 
input (Rm6_1.0) sampled with three different resolu- 
tions. These are the runs Rm6_1.0_12, Rm6_1.0_25 and 
Rm6_1.0_50 listed in Table 2. This is a useful exercise 
because 3d-pdr post-processing requires non-negligible 
time even when run in parallel. Throughout this paper, 
we analyze a coarser resolution than is actually achieved 
by the hydrodynamic simulations. 

Table 3 gives the mean abundance and standard devi- 
ation over all grid points for each of the three sampling 
resolutions. We flnd that differences in the mean abun- 
dances are generally only a few percent and are, with- 
out exception, much smaller than the standard deviation 
of the distributions. The mean gas temperature is also 
fairly insensitive to increasing resolution. 

Figure 3 shows the fractional abundances for a sin- 
gle random sight-line through the cloud. Increasing the 
sampling resolution of 3d-pdr has little effect on the 
calculated cloud chemistry and the abundances of H, II2 
and CO. Different sight-lines exhibit similar good con- 
vergence. The small differences between resolutions im- 
ply that the results should also be similar for simulation 
data with higher base grid resolutions. This comparison 
suggests that in the future it will be possible to follow 
the time-dependent chemical evolution coarsely but accu- 
rately with 3D-PDR. However, for stronger UV flelds, the 
resolution could be more important since the C+/C/CO 
transition will occur further from the boundary. 

3.2. Boundary Convergence 

Some authors have suggested that the details of the 
interior cloud chemistry depend on the resolution of the 
atomic-to-molecular transition. To investigate this issue, 
we compare molecular abundances in the cloud interior 



TABLE 2 
3D-PDR Run Parameters 



Run ID 



Grid (/ X 2563)^ FUV (Gq)^ Final Time (Myr) 



Field Geometry^ 



A/"/ ni..nu (cm-3)« 



Rm6_0.0_12 


1/12 ] 


L 100 


plane-parallel 


12 


200. .10^ 


Rm6_0.5_12 


1/12 ] 


L 100 


plane-parallel 


12 


200. .10^ 


Rm6_1.0_12 


1/12 ] 


L 100 


plane-parallel 


12 


200. .lO'^ 


Rm6_1.0_12i 


1/12 ] 


L 100 


isotropic 


12 


200. .10^ 


Rm6_1.0_12ui 


1/12 1 - 


- h 100 


plane-parallel + isotropic 


12 


200.. 10^ 


Rm6_1.0_25 


1/25 ] 


L 100 


plane-parallel 


12 


200. .10^ 


Rm6_1.0_50 


1/50 ] 


L 100 


plane-parallel 


12 


200. .10^ 


Rm6_1.0_12bf 


1/12 ] 


L 5.7 


plane-parallel 


12 


200. .10^ 


Rm6_1.0_12_48 


1/12 ] 


L 100 


plane-parallel 


48 


50.. 10^ 


Rm6_1.0_12_NC 


1/12 ] 


L 100 


plane-parallel 


48 


0..104 


Rm9_1.0_12 


1/12 ] 


L 100 


plane-parallel 


12 


200. .10^ 


Rm4_1.0_12 


1/12 ] 


L 100 


plane-parallel 


12 


200. .lO'^ 


n300_12f 


1/12 ] 


L 5.7 


plane-parallel 


12 


200. .10^ 



^Input sampling of the simulation data used by 3d-pdr. 
^Magnitude of the UV field in Draines at the box edge. 

^The direction of the field at the boundaries. The field is either a uniform field that is plane-parallel to the box faces, 
isotropic, or a combination of the two. 
^Number of rays. 

^Range of 3d-pdr densities assumed in the calculation. 
^Run uses the same C and O abundances as GIO {xc = 1-41 x lO"'^ and xq = 3.16 x lO"'^). 



TABLE 3 

Mean fractional values at various resolutions 



Snapshot ID^ 



H 



H2 



C+ 



c 



CO 



T(K) 



Rm6_1.0_12 
Rm6_1.0_25 
Rm6_1.0_50 



0.082 (0.16) 
0.099 (0.19) 
0.10 (0.21) 



0.36 (0.20) 
0.35 (0.20) 
0.34 (0.20) 



6.7d-5 (4.0d-5) 
6.7d-5 (4.0d-5) 
6.6d-5 (4.1d-5) 



6.9d-6 (l.Od-5) 
6.8d-6 (l.Od-5) 
6.6d-6 (9.9d-6) 



5.24d-6 (1.5d-5) 45.7 (32.2) 
5.0d-6 (1.5d-5) 45.5 (31.8) 
5.1d-6 (1.5d-5) 44.1 (31.9) 



^Simulation output ID and mass-weighted mean abundances, 
parentheses. 



The standard deviation for each is given in 
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Fig. 3. — Mean fractional H, H2 and CO abundances for a line of 
sight through the cloud center at three different resolutions. The 
resolutions differ by factors of two in the number of grid points. 

for two cloud edge resolutions. Figure 4 shows the same 
sight-line computed with a fixed linear spacing and with 
logarithmically spaced points concentrated at the bound- 
ary. All grid points are assumed to be part of the PDR 
and are treated with the PDR code. We find that the 
abundances in the cloud interior are virtually identical 



despite the very different boundary resolutions. In fact, 
the values computed with coarse resolution vary some- 
what only within one or two coarse cells directly adjacent 
to the boundary. This demonstrates that the chemistry 
in the cloud interior is not sensitive to the edge resolu- 
tion for the densities and FUV field strengths considered 
here and provides further evidence that our lower reso- 
lution 3d- PDR calculations are chemically converged for 
the bulk of the cloud. 

3.3. Ray Convergence 

In order to assess the sensitivity of our results to the 
number of rays, A/^, we compare 3d- PDR calculations 
with 12 (/ = 0) and 48 (/ = 1) rays. In principle, higher 
ray resolution will be more accurate for asymmetric and 
fractal geometries. Figure 5 shows the fractional abun- 
dances for a line of sight through the cloud center. Gener- 
ally, we find good agreement for the two resolutions. The 
H2 and C abundances are almost identical, while some 
differences of up to an order of magnitude are apparent 
for some H and CO points. For H2 and CO, the resolu- 
tion does affect the molecular transition at the boundary, 
where the abundance is lower at higher ray resolution. 
We can understand this by considering the simpler 6- 
ray case for a cell on the domain boundary. Assuming 
that no radiation impinges on the cell from the opposite 
cloud edge, this cell should see 27r sr of the UV field and 
be completely unshielded. However, for 6 perpendicular 
rays, only the ray perpendicular to the boundary will see 
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Fig. 4. — Mean fractional H, H2, C, and CO abundances for 
a line of sight through the cloud center for constant resolution 
(crosses with a solid line) and for logarithmic spacing at the cloud 
edge (diamonds with a dotted line). The right axis indicates the 
gas number density (gray). A uniform 1 Draine radiation field is 
imposed on the cloud from the left {L = pc). 




1 2 roys + ^ 
48 roys O '; 



0.0 
L (pc) 

Fig. 5. — Mean fractional H (solid), H2 (dot-dashed), C (dashed), 
and CO (dottted) abundances for a line of sight through the cloud 
center for Rm6_1.0_12_NC (crosses) with the fiducial ray resolution 
{Afi = 12) and for Rm6_1.0_12_48 (diamonds), which has Afi = 48. 
Here we plot only the abundances for 50 < n < 10^ cm~^, which 
are the 3d-pdr density hmits of Rm6_1.0_12_48. 

the UV field, which results in an angular attenuation of 
47r/6 = 27r/3 sr. Depending on the field strength, this 
may be sufficient to shield the boundary cell from the UV 
field. As more rays are added the angular dependence of 
the field at the boundary becomes better resolved, re- 
ducing the amount of extinction. In Figure 5, we see this 
issue only affects a few cells adjacent to the domain edge 
and does not appear to directly impact the subsequent 
internal cloud chemistry. 

4. RESULTS 
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4.1. Code Comparison: Post-processing vs. In situ 
Calculation 

In this section we compare our results using 3D-PDR to 
the coupled chemical and dynamical method described 
in GIO. There are a few key differences between the two 
approaches. 3D-PDR follows 320 reactions of 33 species 
(including electrons) while GIO follows 218 reactions of 
32 species. We note that these 218 reactions are not an 
exact subset of the 320 followed by 3D-PDR since they in- 
clude more reactions with negative ions. GIO adopts the 
older reaction rates of UMIST99 (Le Teuff et al. 2000) in- 
stead of UMIST07 (Woodah et al. 2007). GIO employs a 
"six-ray" approach (Nelson & Langer 1997, 1999; Glover 
& Mac Low 2007b) to calculate the local attenuated ra- 
diation field whereas 3D-PDR uses A/^ = 12 x 4^ rays 
(in this paper we use 12 rays, i.e. ^=0). Both methods 
include heating due to the photoelectric effect, H2 pho- 
todissociation, UV pumping of H2, H2 formation on dust 
grains, and cosmic ray ionization. However, 3D-PDR also 
includes photo-ionization of neutral Carbon and turbu- 
lent heating. 

Both methods neglect the impact of the gas velocity 
distribution on the chemistry. In practice, the details 
of the velocity field affect the II2 shielding, since the 
H2 photodissociation rate from any given Lyman- Werner 
line is related to the escape probability for that line (see 
Glover & Mac Low (2007a) and discussion therein). GIO 
and previous papers instead adopt a six-ray approxima- 
tion to estimate the shielding, which includes no velocity 
information. 3D-PDR relates the line optical depth to 
an effective linewidth, which is proportional to the root 
mean square of the thermal sound speed and turbulent 
gas velocity. 

In modeling cooling, both methods include emission 
by C, C+ and O fine structure lines, gas-grain collisional 
coohng, cooling by rotational lines of CO, and II2 colli- 
sional dissociation, but 3D-PDR neglects H2 and H2O ro- 
vibrational and OH rotational lines, which are included 
in the one-dimensional code ucl_pdr (Bayet et al. 2010), 
as well as H collisional ionization and Compton cooling. 
However, these lines do not produce significant cooling 
under the conditions considered here, and so neglecting 
them is a good approximation. 

To perform a precise comparison of the two methods 
we apply 3D-PDR directly to the density field of n300. We 
adopt identical C and O abundances {xq = 1.41 x 10~^ 
and xq = 3.16 X 10~^ in all forms relative to hydrogen), 
evolve to the same final time of 5.7 Myr, and apply the 
same stellar FUV field at the boundary. 

Figure 6 shows the H, H2, C, and CO abundances 
binned as a function of effective extinction for the two 
methods. We define the effective extinction following 
GIO: 



A^ 



,eff 




(7) 



where Av,eff is a weighting over the extinctions of all 
A/"^ = 12 rays. For comparison, we include the abun- 
dances for a different simulated density field in order to 
assess the sensitivity to the underlying density distribu- 
tion. In Figures 6 and 7 we consider only densities 200 



cm"^ < n < 10^ cm"^.^ 

We find that the GIO and 3D-PDR results differ the 
most at low extinction (Av,eff < 0.3). These cells are 
near the simulation boundaries, where the impinging ra- 
diation field dissociates the H2. Surprisingly, this tran- 
sition is largely absent in the GIO method, which also 
appears to over-estimate the fraction of atomic hydrogen 
throughout the PDR. At higher extinction, the GIO and 
3D-PDR methods show reasonable confluence between the 
distributions of C and CO. The similarity between the 
n300 and Rm6_1.0_12 distributions illustrates that the 
PDR is not overly sensitive to the underlying density 
distribution. 

We note that the abundance distributions of 
Rm6_0.0_12, which we evolve to 100 Myr, and those of 
Rm6 -0.0-1 2b, which we conclude at 5.7 Myr, are very 
similar. This suggests that these species achieve chemi- 
cal equilibrium by 5.7 Myr (see also Bayet et al. 2009). 
The formation time of molecular hydrogen for gas with 
n = 10^ cm-^ is tform - 10^ n"^ yi c^ 1 Myr (Hollen- 
bach et al. 1971). CO forms rapidly provided A^ > 0.7 
(Bergin et al. 2004) 

Figure 7 illustrates the gas temperature distributions 
in the three cases. GIO reaches a slightly lower tempera- 
ture at low Av, but otherwise the calculations are within 
a standard deviation. The temperature histograms, 
which show the relative number of cells at different tem- 
peratures, are somewhat different. All simulations ex- 
hibit a peak in the temperature distribution at ~ 30 — 50 
K. 

Despite the general congruence of mean properties. 
Figure 8 demonstrates that individual cells may have 
very different fractional abundances. H2 abundance has 
the best point-by-point agreement since it is nearly con- 
stant throughout the domain. The exception occurs at 
the cloud boundaries, where the GIO method does not 
appear to model the PDR regime as well as 3D-PDR. 
The discrepancy is likely due to the lower ray resolution 
of GIO, which causes the H2 fraction to be over-estimated 
(see discussion in section 3.3). GIO mention ray resolu- 
tion as a possible deficit of their 6-ray approach. The 
H abundance shows fairly good correspondence between 
the two methods but does differ by an order of magni- 
tude at some points. The higher H fractional abundance 
at high-Av shown for n300_G10 may be due to turbulent 
mixing, which we do not include in our approach. How- 
ever, since n300 does not include gravity, which would 
reduce turbulent turnover at high-densities, this H frac- 
tional abundance may also be an over-estimate (S. Glover 
private communication). The C and CO abundances pro- 
duced by the two methods are also generally similar with 
the most difference occurring in the range — 5 pc < x < 
pc, which corresponds to gas densities n < 10^ cm~^. 

Since the input densities, grid resolution, and atomic 
abundances are identical, all discrepancies must be due 
to differences in methodology and chemical assumptions. 
Although the magnitude of the abundance variation ap- 
pears quite large, such differences are consistent with 
those typically found between PDR codes (Rollig et al. 

^ We find that Figures 6 and 7 appear very similar assum- 
ing a lower cutoff of n > 50 cm~^ (e.g., Rm6_1.0_12_48 and 
Rm6_1.0_12_NC). The main result of including lower densities in 
the PDR calculation is that the low-Av gas becomes increasingly 
(and inaccurately) warm. 



2007). 

4.2. Dependence on Dimensionality: ID versus 3D 

Figure 9 shows the extinction distribution as a function 
of the local gas number density. As shown by GIO (e.g., 
their Figure 14), gas extinction and density are only very 
weakly correlated. Overall, extinction is more strongly 
correlated with position within the cloud than with den- 
sity. Figure 9 exhibits two populations of points: those in 
the interior and those spatially within < 1% the bound- 
ary, which have distinctly lower extinction and stand out 
in the regime 3.5 < logn < 4. These boundary cells ap- 
pear to be missing from the results of GIO. We discuss 
them further in section 4.3. 

Figure 10 illustrates how the H2, C+, C, and CO frac- 
tional abundances depend upon extinction, UV field, and 
location within the cloud. Both H2 and CO, which are 
the most sensitive to the extinction, appear to behave 
differently at the cloud edge and in the interior. Some 
of the discontinuity in the distribution is likely artificial 
since better edge resolution, as shown in Figure 3, would 
join the two populations more smoothly. 

The left column plots verify that the the local extinc- 
tion and UV field magnitude are completely correlated 
(i.e., the extinction indicated by the color-scale varies 
completely linearly with the magnitude of the UV field). 
Relative to extinction, proximity to the cloud edge has 
a weaker influence since the abundances depend upon 
the column density along each ray, which varies as the 
density distribution. 

Figure 11 illustrates how abundances correlate with 
the gas temperature. Temperature varies smoothly with 
both UV and Av,eff- For C, a discrete region of boundary 
cells becomes apparent, which was previously degenerate 
with the cells near but not abutting the boundary. These 
boundary cells show up as a slight offset in temperature 
for Av,eff < -0.5. 

The shape of the abundance-UV field distributions de- 
pends on the range of underlying densities in the PDR. 
Figure 12 shows the abundance distributions over plot- 
ted with lines showing the abundances computed for sim- 
ple ID models. The ID models assume constant density 
along the line-of-sight and an incident 1 Draine UV field 
at one end. These curves illustrate that the range in 
abundance for any given UV field simply depends on the 
range in local gas density for a given UV field. Since 
we define the PDR region as those points with densi- 
ties n = 200 cm~^ to 10^ cm~^ the curves with these 
densities correlate well with the data of the 3D simula- 
tions. Thus, while turbulence dictates the distribution 
of densities and hence the fraction of cells within a given 
density range, the abundance distribution is set by chem- 
istry and the details of the species response to the local 
UV field. In summary, although the determination of 
the extinction at a particular point within the volume is 
a three-dimensional problem, we find that once the lo- 
cal UV field is computed, the resulting abundances are 
nearly identical to those derived from a ID model. 

4.3. Dependence on Physical Parameters 

In this section we investigate the sensitivity of the 
chemistry to the bulk simulation properties. In a self- 
consistent treatment of molecule formation, the distri- 
bution of shock properties (e.g. the post-shock densities 
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Fig. 6. — Mean fractional H, H2, C, and CO abundances as a function of log extinction, where the error bars indicate the standard 
deviation in each bin. The triangles show the results of the n300 simulation as performed by S. Glover. The diamonds indicate the results 
for 3D-PDR assuming the same n300 density distribution evolved to the same time (5.7 Myr) using the same xc and xq abundances. 
The squares indicate the results for Rm6_1.0_12b, which has a different underlying density distribution. In all cases, only gas with 200 
cm~^ >n< 10^ cm~^ is included in the averages. 
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Fig. 7. — Mass- weighted temperature distribution versus extinction (left) and temperature distributions (right) for the same runs shown 
in Figure 6. 
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Fig. 8. — Fractional abundances computed for a line-of-sight 
through the cloud center of simulation n300 for the GIO method 
(narrow lines) and 3d-pdr (thick lines). The gas number density 
(right axis) is indicated by a gray, solid line. For comparison, no 
cutoff is applied to the 3d-pdr calculation. 
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Fig. 9. — Log extinction versus number density for Rm6_1.0_12. 
The colorbar indicates the minimum distance to a cloud boundary 
in pc. 

and temperatures) could imprint an observable signature 
in the measured abundances. For example, molecules 
such as CH2, HCO+ and OH are directly sensitive to 
turbulent density fluctuations Kumar & Fisher (2013) 
and would likely vary as a function of simulation Mach 
number. In contrast, the abundance of H2 and CO pre- 
dominantly depends on the amount of shielding from the 
UV radiation field (Bergin et al. 2004). A parcel of gas 
embedded within a completely smooth {Ay > 0.7) cloud 
will be well-shielded from the UV field, whereas a parcel 
of gas in a highly fractal cloud will have a high probability 
of having a sight-line with low extinction. Consequently, 
we can expect that the morphological distribution of the 
gas will have some effect on these abundances. 

Figure 13 shows the mass- weighted abundances for 
simulations with two different Mach numbers at various 
evolutionary times. All ORION simulations have the same 
mean density such that apparent differences are directly 
due to variations in the gas morphology. 



Due to the relatively high simulation mean-density, we 
find that the H2 abundance is fairly insensitive to changes 
in the density distribution caused by gravity. Glover & 
Mac Low (2007a) found that the H2 fraction increased 
with gravitational collapse for a smooth density distri- 
bution. However, in their case the initial mean densi- 
ties were much lower than the initial density of our runs. 
Thus, in our simulations the gas is predominantly molec- 
ular at all times since most gas parcels are well-shielded 
with or without gravity. 

We find that C and C+ abundances vary by less than 
20%. The CO abundance demonstrates the largest vari- 
ation and declines by more than a factor of two as the 
gas becomes self-gravitating. Much of this effect is due 
to a non-negligible fraction of the mass becoming con- 
centrated in small, dense and well-shielded volumes that 
have fixed, maximal 10 ~^ abundance, which decreases 
the extinction in the remainder of the volume (see Fig- 
ure 2). Differences of '^30% in 3D Mach number have a 
relatively small impact (less than a factor of 2) on the 
total mean abundance but CO does show some sensitiv- 
ity with abundance increasing with Mach number. H2 
and C+ fractions decrease slightly with increasing Mach 
number. 

For reference. Figure 13 includes the mean abundances 
of n300. In order to compare the PDR results we only 
consider n300 gas with densities exceeding 200 cm~^.^ 
In this case, the abundance differences are dominated by 
the different mean extinctions. n300 has Av ~ 0.04 while 
Rm6 has Av ~ 0.02, which results in slightly lower mean 
C and CO abundances. The n300 mean H2 abundance 
is lower than for those computed for Rm6; however. Fig- 
ure 4 in GIO exhibits a higher H2 abundance (~ 0.98) 
for a simulation with the same mean density but lower 
magnetic field. This suggests that there is a morpholog- 
ical component to the H2 abundance difference that is 
related to the magnetic field strength (S. Glover private 
communication) . 

Figure 14 illustrates the CO distribution as a function 
of gas temperature. As gravity infiuences the gas dis- 
tribution, the number of cells with high CO abundance 
and cold temperature (10 < T < 20) increases. This 
is related to the volume filling factor of the dense gas, 
which decreases as gas becomes more concentrated in 
dense, collapsing regions. The shape of the temperature- 
abundance distribution is otherwise roughly constant 
with Mach number and time. 

In all panels, points near the edge of the simulation 
box comprise a distinct swath of high-temperature/low- 
abundance points. We color points within 2% of the 
edge red to highlight this dichotomy. This region directly 
corresponds to the low- Ay, mostly atomic region at the 
boundary. 

4.4. Dependence on External Radiation Field 

In order to investigate how abundance depends on the 
external radiation field, we consider two 3D-PDR calcu- 
lations with external fields each with a magnitude of 
1 Draine but with different vectors. Run Rm6_1.0_12i 
has an external isotropic radiation field, while Run 

^ The n300 mean abundances are fairly similar whether or not 
the lower density gas is included because the averages are mass- 
weighted. 
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Fig. 10. — Left: Rm6_1.0_12 fractional abundances as a function of log UV radiation field where the colorbar indicates the extinction. 
Right: Rni6_1.0_12 fractional abundances as a function of log extinction (right), where the colorbar indicates the minimuni distance in pc 
to the cloud boundary. 
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Fig. 12. — Log abundance versus log UV field for run Rni6_1.0_12, where the color scale indicates the log of the mean density for a given 
abundance and local field. The gray solid lines show the values from ID runs with constant density and incident 1 Draine UV field. The 
lines increase in density for line thickness going from thin to thick with values of 200, 500,10"^, 2 x 10*^, 5 x 10"^, 10^ cm""^. 
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Fig. 13. Mean H, H2, C+, C, CO abundances for Rm6_0.0_12, Rm6_0.5_12, Rm6_1.0_12 (red), Rm9_1.0_12 (purple), Rm4_1.0_12 (green) 
and n300 using the GIO method (blue). On the right abundances have been normalized to the values of Rm6_1.0_12. 
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Fig. 14. — CO abundance as a function of gas temperature for different times and Mach numbers. The cells that are within 2% of the 
box boundary are colored red. 
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Fig. 15. — Schematic of the incident radiation field used in 
the 3D-PDR runs. The cube represents the entire computational 
domain. The thin solid lines represent the boundaries of the 12- 
ray healpix structure as they are emanated from a point placed 
in the center of the computational domain. The direction of the 
isotropic radiation field is opposite to the direction of each healpix 
ray as shown in the figure. The additional dashed lines on the left 
represent the direction of the plane-parallel radiation field added 
in run Rm6_1.0_12ui. 

Rm6_1.0_12ui has an external field that is a superposition 
of a half Draine isotropic field and a half Draine uniform 
field (i.e., a field that is plane parallel to the simulation 
boundary at all faces). Figure 15 illustrates the incident 
radiation field geometry for the two cases. 

Figure 16 illustrates that by simply changing the field 
incidence the internal point-by-point UV distribution is 
very different. The figure includes only points with a 
net field greater than 0.5 Draines; these points are ones 
which feel both the isotropic and plane-parallel compo- 
nents of the incident field and thus display the maximum 
difference. Figure 17 shows the effect of the field differ- 
ences on the H, H2, C, and CO fractional abundances. 
Since CO abundance depends mainly on the local UV 
field and these distributions are distinct, it is unsurpris- 
ing that individual abundances change by as much as 
50% for different field configurations. Likewise, C and 
H are strongly affected by the field distribution. Figure 
17 shows that the mixed field simulation has fewer high 
UV points and more lower UV points, which is consis- 
tent with the elevated C and CO abundances displayed 
in Figure 17. Since the molecular hydrogen abundance 
is nearly constant within the cloud, the field configura- 
tion at the boundary has little effect. The higher density 
gas (n > 10^ cm~^), which is well self-shielded by defi- 
nition, is also largely insensitive to field changes of this 
magnitude. 

In summary, even a modest change in the UV field in- 
cidence reinforces the conclusion that three-dimensional 
PDR treatment is preferable to a one-dimensional treat- 
ment for complex or non-symmetric problems. We ex- 
pect differences to be more significant for larger exter- 
nal field variations and for the inclusion of internal UV 
sources, i.e., protostars. 

5. CONCLUSIONS 

We use 3D-PDR in combination with hydrodynamic 
molecular cloud simulations to explore the importance 
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Fig. 16. Distribution of UV field for the Rm6_1.0_12i run which 
has an isotropic-only external 1 Draine field (blue, dotted) and the 
Rni6_1.0_12ui run which has a 0.5 Draine isotropic and 0.5 Draine 
uniform external field (solid, purple). Only the points that have 
a net field greater than 0.5 Draines are plotted. The two 3d-pdr 
runs use the same input density distribution. 

of dimensionality in PDR chemistry, to consider com- 
plex gas morphologies and to compare with prior results 
using an in situ astrochemistry treatment. 

First, we demonstrate that our results are robust as a 
function of grid sampling and edge resolution. In fact, we 
find that the interior cloud abundances are remarkably 
insensitive to the resolution of the atomic to molecular 
transition at the cloud boundary. 

We obtain reasonable agreement between the GIO in 
situ and 3D-PDR approaches for the C and CO abundance 
distributions. This is because molecules such as CO and 
H2 are not particularly sensitive to the dynamical his- 
tory of the gas but instead depend predominantly on the 
local radiation field. The two approaches differ the most 
for H and H2 abundances near the cloud boundary and 
for cells that have low-extinction. For example, in GIO 
hydrogen is either entirely molecular or fully dissociated. 
This discrepancy appears to result from differences in 
the methodologies rather than chemical details, and we 
assert that the treatment of 3D-PDR should be more ac- 
curate in transition regions. 

We demonstrate that morphological differences due 
to cloud Mach number and evolutionary time can pro- 
duce significant differences in the abundance distribu- 
tions. While this may be difficult to observe directly 
since point by point abundances are difficult to infer, 
it may indirectly impact the properties of the observed 
molecular emission lines emerging from the cloud. 

Finally, we find that a relatively modest change in the 
external UV radiation field produces large changes in the 
chemical abundances. This supports the finding by B12 
that three-dimensional treatment is crucial for complex 
and non-symmetric problems. 

In paper II, we plan to implement several improve- 
ments to our method. First, we will relax our simple 
abundance approximations at densities > 10^ cm~^ and 
instead couple 3D-PDR to a one-zone gas-grain chemi- 
cal network code. This will allow us to include molec- 
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ular freezeout onto dust grains, gas-grain thermal cou- 
pling, and a more extensive chemical network. Second, 
we will investigate the effect of embedded UV sources on 
the chemical distribution. The ORION simulations that 
we analyze here contain detailed information about the 
masses and accretion rates of embedded protostars that 
we have neglected in this study. In addition to these 
changes, work is ongoing to couple 3D-PDR to the pho- 
toionization and radiative transfer code MOCASSIN. 



10000 



1000 

n (cm"^) 



10000 



and CO abundances for Rni6_1.0_12i and Rni6_1.0_12ui (the same runs 
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APPENDIX 

HEATING RATES 

There are four main contributions to the local heating at each domain point. First, there is photoelectric heating, 
which is produced by UV photon interactions with dust grains and PAHs, and which dominates near the cloud surface. 
Second, there is cosmic-ray ionization heating, which is set by the standard cosmic-ray density and becomes dominant 
deeper into the cloud. Third, there is chemical heating as the result of various exothermic reactions. Finally, there 
is turbulent heating, which is due to energy dissipation through shocks and depends on the turbulent outer scale 
(e.g., cloud size) and turbulent Mach number. Figure Al shows the distribution of heating rates for each of these 
contributions. Photoelectric heating dominates in most cases, and the turbulent heating generally provides the smallest 
contribution. If the turbulent heating is proportional to v^jj^^^/L, then for the 2pc simulations here it should range 
from '^turb/^ = 1.5 x 10~^ — 1.4 x 10~^ cm^ s~^, which brackets the constant value, 6.5 x 10~^ cm^ s~^, we adopt in 
our calculations. We direct the reader to Pan & Padoan (2009) and Kumar & Fisher (2013) for additional discussion 
and modeling of heating due to turbulent dissipation, intermittency, and shear fiows. 
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Fig. A1. — Distribution of heating rates in Rni6_1.0_12 as a function of effective extinction for the total photoelectric (blue, top contours), 
cosmic-ray (green, top-middle contours), chemical (black, bottom-middle contours) and turbulent heating (red, bottom contours) in each 
cell. The contours correspond to the number of cells, n, with a given extinction and heating rate: n =100 (inner contour), n = 30 (middle 
contour), and n = 10 (outer contour). 
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